if (!require("pacman")) install.packages("pacman")
pacman::p_load(dplyr, tidyr, ggplot2)

# Load and filter data

dat <- readRDS("data/zij.rds") %>%
  filter(yob > 1959 & yob <= 1975)

# Define outcomes

outcome_vars <- c("DDR_kapitalismus", "vereinigung_dagegen", "DDR_reformsozial")

# Reshape and recode data for plotting

plot_df <- dat %>%
  pivot_longer(
    cols = outcome_vars,
    names_to = "outcome",
    values_to = "dv"
  ) %>%
  mutate(
    outcome = recode(outcome,
      "DDR_kapitalismus" = "Support capitalism",
      "vereinigung_dagegen" = "Oppose reunification",
      "DDR_reformsozial" = "Support reformed socialism"
    ),
    yob = as.numeric(yob)
  )

# Calculate yearly means by group

plot_df_sum <- plot_df %>%
  group_by(gender, post_gender, outcome, yob) %>%
  summarise(m = mean(dv, na.rm = T)) %>%
  ungroup()

# Generate plot

p <- ggplot(plot_df, aes(x = yob, y = dv, group = post_gender, col = gender)) +
  geom_smooth(method = "lm", se = FALSE) +
  geom_point(data = plot_df_sum, aes(x = yob, y = m, color = gender)) +
  geom_vline(xintercept = 1970.5, linetype = "dotted") +
  scale_color_manual(values = c("black", "grey")) +
  scale_x_continuous(breaks = 1960:1975) +
  facet_wrap(~outcome, nrow = 3, scales = "free_y") +
  theme_bw() +
  theme(
    legend.position = "bottom",
    axis.text.x = element_text(angle = 45, hjust = 1),
    axis.title.x = element_blank()
  ) +
  labs(col = "", y = "Mean of DV")

p
